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In this research, a simplified uniaxial strain -controlled creep damage law is 
deduced with the use of experimental observation from a more complex 
strain-dependent law. This creep damage law correlates the creep damage, 
which is interpreted as the density variation in the material, directly with 
the accumulated creep strain. Based on the deduced uniaxial strain-controlled 
creep damage law, a continuum mechanical creep rupture analysis is carried out 
for a beam resting on a high temperature elastic (Winkler) foundation. The 
analysis includes the determination of the nondimensional time for initial 
rupture, the propagation of the rupture front with the associated thinning of 
the beam, and the influence of creep damage on the deflection of the beam. 
Creep damage starts accumulating in the beam as soon as the load is applied, 
and a creep rupture front develops at and propagates from the point at which 
the creep damage first reaches its critical value. By introducing a series of 
fundamental assumptions within the framework of technical Euler -Bernoull i type 
beam theory, a governing set of integro-dif ferential equations is derived in 
terms of the nondimensional bending moment and the deflection. These 
governing equations are subjected to a set of interface conditions at the 
propagating rupture front. A numerical technique is developed to solve the 
governing equations together with the interface equations, and the computed 
results are presented and discussed in detail. 
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1. INTRODUCTION 


Tertiary creep involves the process of fracture leading ultimately to 
complete failure, and is associated with local reduction in cross-sectional 
area and more importantly with the nucleation and growth of voids and 

microracks along grain boundaries. This failure mode leads to eventual 
collapse of a structural component and is known in the literature as creep 
rupture or stress rupture. In order to meet the demands of designers and 

engineers concerned with the safety of equipment operating at elevated 
temperatures, researchers in recent decades have conducted extensive creep 
rupture experiments from which they hope to extract some useful 
"extrapolation" parameters. Such parameters are inevitably limited by the 

laboratory-allowed time scale and by the usual scatter of the empirical data, 
but they are employed to estimate the appropriate stress and temperature 
requirements for the practical service lives of equipment in operation. 
Amongst such extrapolation parameters methods are the ones of Larson-Miller 
(1952) [1], Manson-Haferd (1953) [2], Orr-Sherby-Dorn (1954) [3], and many 

others. Manson and Insign [4] have presented an interesting review on the 

progress in extrapolation procedures for creep rupture; an excellent 
discussion of these is also given in the text by Conway [5]. 

In parallel with the development cited above, other researchers including 
some metallurgists have attempted to define and quantify a suitable variable 
which describes the damage state and measures the extent of damage in 
materials undergoing creep. The major hurdle in this line of research is the 
manner by which ones bridges the gap between the scalar damage variable 
obtained by macroscopic creep testing and the microscopic processes involved 
in the nucleation and growth of voids and microcracks at grain boundaries. 
Such variables are expected to be able to characterize the damage state from 
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the physical and quantitative points of view, and furthermore to provide a 
useful tool for analytical modelling via continuum mechanics. Amongst such 
approaches are Robinson's linear cumulative creep damage law (1952) [6]* 
Hoff's ductile creep rupture theory (1953) [7], Kachanov's brittle rupture 
theory (1961) [8], Robotnov's coupled damage creep theory (1969) [9], and many 
other modified theories such as the one due to Leckie and Hayhurst (1974) 
[10], Comparative studies of the various theories may be found in [11-13]. 
Recently, scientists have observed a close relation between density change and 
the nucleation and growth of voids and microcracks associated with creep 
damage in polycrystalline materials. Extensive efforts have thus been made to 
identify and quantify creep damage in terms of the density variation which is 
attributed to cavitation in a creeping material. Following this concept, 
Piatti et al [14] developed a refined experimental technique to measure the 
density variation for use as a definition of creep damage. Using data 
obtained in this manner for steel, Belloni et al [15,16] proposed a 
statistically-based damage law in a complicated power law form similar to the 
one presented in Woodford's parametric study of creep damage [17]. 

Because of its inherent mathematical complexity, the creep damage law 
proposed in [15,16] is somewhat inconvenient for analytical treatment within 
the framework of continum creep damage mechanics. In addition, some 
arbitrariness remains in the determination of the material constants appearing 
in this damage law (see [18]). Accordingly, the first task in this work is to 
obtain a simplified yet still useful damage law. This task is addressed in 
Section 2 where we argue first from the microscopic point of view that density 
variation certainly Is a proper Index of damage In a material undergoing creep 
deformation. We then propose a simplified uniaxial strain-controlled damage 
law by introducing some assumptions based on experimental observation 
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associated with the original damage law, and this strain-control led damagelaw 
is demonstrated to be closely related not only to the original damage law but 
also to Kachanov's damage law (see [8]). We conclude Section 2 with the 
observation that, whereas a typical boundary value problem suffices to 
represent the problem in “the first stage of creep damage", we encounter in 
“the second stage of propagation of the rupture front" a moving boundary 
problem similar to the Stefan problem in heat conduction [19]. 

Utilizing the above strain-controlled creep damage theory, we present in 
Section 3 a continuum mechanics model for the creep rupture analysis of a beam 
resting on a high temperature elastic Winkler foundation which generates a 
prescribed thermal gradient in the thickness direction. Based on technical 
Euler-Bernoulli -type beam theory, we derive in Section 3 a set of governing 
differential equations for a region with a moving boundary (rupture front) 
which is prescribed by a set of interface equations. Owing to the inherent 
nonlinearity of the problem, closed form solutions generally do not exist. 
Accordingly, a successful treatment of the problem requires the application of 
a suitable numerical technique which is then presented in Section 4. In the 
latter part of Section 4, we consider a simple case for which a closed form 
solution does exist. We then present detailed numerical results for the 
problems in which temperature gradient is taken into account and the 
foundation is either present or absent. The results consist of the 
nondimensional forms for bending moment, deflection, and the geometric shapes 
of the rupture front. 


2. STRAIN-CONTROLLED CREEP DAMAGE 

2.1 Creep Damage Law Under Uniaxial Stress 

Virtually all load-bearing structural components operating at elevated 
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temperatures undergo the typical 3-stage creep phenomenon. Various 
phenomenological Interpretations of the creep process have been devised, 
usually employing the concept that creep Is essentially a competition between 
strain-hardening and recovery [20]. It Is well understood that at elevated 
temperatures a crystalline solid may deform In accordance with several 
mechanisms such as dislocation creep and diffusion creep. Each such mechanism 
Is most active Is some range of stress and temperature [21], such that within 
certain regions of the stress-temperature space one mechanism Is said to 
dominate the others. The pictorial maps constructed by this concept are known 
as Ashby's deformation-mechanism maps [21,22]. Raj and Ashby [23] have 
pointed out that the creep mechanisms mentioned above are In fact an 
"accommodation process" for grain boundary sliding. When a shear stress 
causes sliding to occur at a generally nonplanar grain boundary, some 
accommodation process (such as diffuslonal flow or plastic flow) Is necessary 
to heal the crystalline structure at the deviation of the boundary from a 
perfect plane. In the event that this accommodation process does not develop 
fully at a boundary deviation during sliding, an "incompatibility" results In 
the form of voids and wedge cracks along the grain boundaries which are 
oriented roughly perpendicular to the tensile axis. As the material Is 
strained further the coalescence of voids and cracks eventually leads to 
intergranular creep fracture. Clearly, as the cavity volume Increases during 
the process of tertiary creep and eventual fracture, the material dilates. In 
this section, we shall focus on a strain-controlled constitutive contlnum 
damage law based on this close relation between creep damage and cavitation 
Induced dilation In materials. 

The type of damage described above is associated with power-law or 
dislocation creep [23,24], Steady dislocation creep under constant uniaxial 
tensile stress a 0 is found experimentally to obey the constitutive 
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relation[25] 


c = A(T) o" (1) 

s ° 

in which n is the constant stress power. The reciprocal viscosity coefficient 
A(T) is expressed as the Arrhenius equation 

A(T) = A* exp( -AH/RT) (2) 

★ 

where A is the relatively temperature insensitive pre-exponential cofficient, 
AH the activation energy for creep, R the gas constant and T the absolute 
temperature. Equation (1), which is also known as Norton's steady creep law, 
will be employed to describe the creep deformation process in the problem 
considered later in this work. 

From the phenomenological point of view, creep rupture can be separated 
into two categories. Failure at high stress and low temperature Is 
characterized by pronounced lateral contractions and the first continuum model 
for this process is known as Hoff's ductile creep rupture theory [7]. On the 
other hand, low stress levels together with high temperatures result in 
brittle type of rupture with little lateral contraction, and the first 
phenomenological theory for this process was formulated by Kachanov [8]. We 
shall not consider Hoff's theory further here (ample discussion is given in 
[13,26]), but we shall now review Kachanov's theory briefly. Kachanov defined 
the damage variable « for a one-dimensional test specimen in accordance with 

A -A A 

° e e 
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where A 0 and A fi are, respectively, the original and the effective cross 
sectional areas carrying the load. Clearly, cavitation creates new Internal 
surface area which In turn reduces the effective cross-sectional area carrying 
the load. Thus, the material In Its virgin state has the damage u equal to 
zero, while the damage in a completely deteriorated material approaches 
unity. A power law for the damage-rate was postulated by Kachanov for 
variable one-dimensional stress as 

4 - fe] v (3) 


where C,v are material constants, and where C may be temperature dependent. 
Assuming that the material is initially undamaged, integration of the above 
equation gives 


1 - (l-w) v+1 = C(Uv)jJ o v (t')dt' (4) 

As pointed out earlier In this section, a close connection exists between 
creep damage and the cavitation Induced dilation of a material. Belloni et al 
[15,16] have employed material density variation as the measure of damage in a 
creep material using refined techniques. They proposed a damage law at 
constant stress o 0 in the power form for the uniaxial tension test 

a y 6 

D = c c a t (5) 

C o 
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where D= -A p/p 0 , p 9 is the density of the material in the virgin state, and 
A p is the change in density due to the volume dilation of the material. In 
eqn.(5) e c denotes the creep strain, and c, <*,y,A appear to be relatively 

insensitive to temperature, but c is highly temperature sensitive. In analogy 
with Kachanov's damage variable u , the damage D has value equal to zero in 
the virgin state and is equal to a critical value at rupture D r , which is a 
material constant. Employing statistical regression techniques, Belloni et al 

were able to correlate their experimental data with damage law (5). A close 

Inspection in-eqn.(5), however, reveals that some arbitrariness exists in the 
determination of the material constants (for details see [18]). This 
arbitrariness is a consequence of treating c and a 9 as Independent state 

variables In eqn.(5), without considering the constitutive creep law. One 
possible way of eliminating this arbitrariness is outlined by the sequence of 
simplifications given below. 

First, based on the findings in [16] and related work [26,27] we shall 
make the simplifylg assumption 


Y = An (6) 

It will be shown later that eqn.(6) together with v=n in eqn.(4) establishes 
an equivalence between Kachanov's formulation and the current one. We further 
assume that steady creep as described by eqn.(1) completely dominates the 
deformation behavior, i.e., the material is non-Newtonean viscous. A 
combination of eqns.(l) and (5), together with assumption (6), then gives 


c 

D = 

— r 
A(T) 


a+A 


(« t) 


s 


( 7 ) 
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or 


c <*+6 

c 

X s 

A(T) 


( 8 ) 


in which c s is the creep strain under the steady creep condition. It has 
been found [28] that, for the rupture mechanism considered here, the product 

of steady-creep-rate, c , and the time-to-rupture, t D , is a constant, i.e. 

S K 


c t 
s R 


= C 


MG 


where C u „ is known as Monkman-Grant constant which has the dimension of 
Mb 

strain. This relationship holds true for a wide range of temperature and 
stress. Therefore, at rupture eqn. (7) gives 


D = 


r 


X 

A(T) 


(C ) 
MG 


o+8 


(9) 


where D r is the critical value of damage at rupture. Belloni et al's data 

[15] showed that the critical value of damage in the high temperature range is 

relatively insensitive to temperature. The temperature independent character 

of both C ur and immplies that the B/A(T) 4 in eqn. (9) must also be 
Mb r 

temperature independent. Thus 
c 

= c 

Z * 

A(T) 


where c 0 is a temperature independent material constant. Substituting this 
back into damage law (8) we then get the simplified form 


a+6 

0 = C c 
° S 


(10) 
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Although a significant simplification has been obtained, damage law (10) 
is still physically plausible. Note that, although damage is an explicit 
function of strain along, it is an implicit function of temperature and stress 
via creep constitutive law (1). In accordance with eqn.(10), a material 
exposed to stress experiences damage directly related to the creep strain, and 
rupture occurs as the available creep ductility is exhausted. Hanna and 
Greenwood [29] showed, for copper with pre-nucleated cavities subjected to low 
stress and with the creep rate linearly related to the stress, that 


- Ot G (11) 

P C 

o 

Although a surprising analogy appears to exist between eqns. (10) and (11), 
conclusions may not be easily drawn on the material constants in eqn.(10). 
However, it does appear very reasonable to postulate that creep damage, as 
measured by density variation, be expressed explicity as a function of creep 
strain. 

In many engineering practices, however, the stress may be varying with 
time due to effect such as load variation and stress relaxation. The 
extension of the original creep damage law, eqn.(5) , to the case of time 
dependent uniaxial stress has been presented in [30]. In the case of 

simplified eqn.(10) it suffices to employ the integral form of creep strain 
for variable stress, and thus integrating e $ = A(T)a(t) n we obtain 


D(t) 



o n (t')dt 



( 12 ) 


Here, the reciprocal viscosity function, A(T), Is retained Inside the Integral 
sign, in order to allow for the situation in which the temperature varies with 
time. If v=n In Kachanov's theory [see eqn.(4)], eqns. (4) and (12) assume a 
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very similar form; a more detailed comparison of strain-dependent theory with 
Kachanov's approach Is given In [31]. 

2.2 Propagation of a Creep Rupture Front — The Moving Boundary Problem 

A nonuniform state of stress may be Introduced by the Irregular geometry 
of a structure, nonhomogeneous material properties, and nonuniform external 
loads. Under such circumstances the creep damage within the structure would 
be a function of the space coordinates in addition to time. Creep damage 
starts accumulating In the structure as soon as the loads are applied. As 
time elapses, the creep damage at some point within or on the surface of the 
structure would first reach the critical value, D r , at which rupture takes 
place. This initial rupture time, t. , is determined in accordance with 

eqn.(12) as . 

, *1 “ +s 

D r * c „ { J A(T) <r n (t')dt'} (13) 

o 

A rupture front then develops generally as a smooth surface, and starts 
propagating through the structure until the entire structure collapses at some 
time tjj. It is readily seen that the lifetime of a structure may be divided 
into two time intervals or stages, l.e., 0 < t < tj and tj< t < tjj. In 
the first stage 0 < t < tj, which has been termed the stage of latent 
failure by Kachanov [8] or the Incubation period by Johnson [32], the creep 
damage is assumed to be less than the critical value (D^) everywhere in the 
structure. In the second stage tj < t < t^, which has been termed the 
stage of propagation of rupture, a rupture front l along which 

D = 0 r , (14) 

travels through the structure and complete collapse occurs at t^. 
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A condition on the direction of travel for the rupture front £ may be 
obtained by taking the total time derivative of eqn.(14). Accordingly, we 
obtain 

3D ,3D dN 

+ as 0 

at aN dt 


in which N designates the coordinate normal to the rupture front. Similarly, 
the geometry of the rupture front l is constrained by 


3D 3D 1 

3t dxj dt 


0 


(15) 


in which the are the space coordinates. 

This type of problem, more generally called moving (free) boundary 
problem, is well-known In heat conduction (19] with phase change and is known 
as the Stefan problem. Although the Stefan problem and the creep rupture 
problem share analogous mathematical characteristics, there are some 
significant physical differences. Instead of the temperature profile of the 
Stefan problem, we are now more concerned about the mechanical behavior of the 
structure, such as the coupling between the stress redistribution and the 
speed of the moving boundary (rupture front). Owing to the inherent 
nonlinearity of the problem, closed form solutions generally do not exist for 
moving boundary problems with a finite damain. Accordingly, a successful 
treatment of such problems will require the application of a suitable 
numerical technique. A thorough discussion and comparison of numerical 
methods currently used for moving boundary problems is given in [33]. Details 
of the numerical technique which we shall choose will be disclosed in 
subsequent sections as the need arises. 
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3. CREEP RUPTURE IN A BEAM ON A WINKLER FOUNDATION 


3.1 Statement of the Problem 

The beam problem to be studied Is depicted in Fig. la. We consider a 
beam continuously supported by an elastic Winkler foundation, which exerts a 
restoring force as the beam deflects under the action of a distributed lateral 
load. Since the foundation is at an elevated temperature, a prescribed 

thermal gradient is assumed to exist in the z-dlrection (thickness) of the 
beam. It Is assumed that this prescribed temperature distribution through the 
thickness of the beam is Independent of time during the deformation and 
rupture processes. The physical model used to analyze the problem Is shown In 
Fig. lb. Here, the elastic foundation is modelled as an infinite series of 
infinitesimal springs with an elastic constant K [34], i.e. as an elastic 
Winkler foundation. Creep deformation starts to accumulate in the beam as 
soon as the lateral load Is applied. In geophysical research this type of 
flexure model has recently yielded some Interesting results on lithospheric 

flexures (eg. see McMullen et al [35]), where the temperature variation with z 
is due to the geothermal gradient and the Winkler foundation is due to the 
underlying mantle. In addition to the creep deformation, we shall also 
consider the effects of creep damage using the concepts previously developed. 
In brief, it is our major goal here to explore the propagation of a creep 

rupture front in a non-isothermal beam under distributed lateral load. During 

the second stage of damage the beam is thinning in a non-uniform manner, and 
accordingly the cross-section of the beam Is not constant (see Fig. lc). It 
will be seen later that a moving boundary problem is encountered as a 
consequence of this thinning behavior. 

The problem presented here is extremely complex in nature. In order to 
reduce the mathematical difficulties somewhat, we present below a series of 
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simplifying assumptions. Firstly, we assume that the material in the beam 
obeys the Norton law of steady creep, with viscosity dependent on the 
prescribed temperature gradient. Although the beam is of non-uniform 
cross-section during the second stage of damage, we assume that technical 
Euler-Bernoul li-type beam theory Is valid throughout the entire process of 
creep damage. We also restrict our consideration to the case of small 
deformations and small rotations. Furthermore, we assume that no major cracks 
form In the unruptured segment of the beam during the process of rupture, and 
thus the effects of stress concentration at crack tips are excluded from the 
current study. Finally, we assume that the shear stresses are negligibly 
small when compared with the axial stresses due to flexure. 

3.2 Mathematical Formulation of the Problem 

The constitutive law governing the creep deformation in the beam Is 
assumed to be of the Norton type [eqn. (1)]: 

c c = A(z)o n (16) 

Here the stress state a may vary with time as well as with the x- and 
z-coordinates. Note also that the reciprocal viscosity coefficient function 
A(z) is an implicit function of z via the temperature distribution (see Fig. 
lc). 

The geometry of the beam is shown is Fig. lc; it has a rectangular 
cross-section of width b and thickness h and the length of the beam is 2L. 
For simplicity we will consider symmetric loading and thus only symmetric 
deformation in this work, and therefore only half of the span of the beam need 
be considered. Employing Euler-Bernoul li-type beam theory with h constant, we 
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may derive the expression for stress in terms of the bending moment M as 


a 



(17) 


where e c Is the distance to the neutral axis (marked as N.A. in Fig.lc). 
Also, the governing equation in the bending moment M is obtained as 


. K( L)H = (18) 

8x«0t ,6 «X 4 3t 

o 

where P Is the applied lateral load, and we have Introduced the notation for 
flexural rigidity 


So 



(19) 


The R.H.S. of eqn.(18) vanishes if we assume that the applied lateral load 
P(x,t) is expressed mathematically in the form 


P(x,t) = P 0 f(x) H(t) (20) 

where P 0 is the maximum load at x=0, f(x) is the symmetric shape function and 
H(t) represents the Heaviside unit step function. For a viscous material 
governed by eqn. (16) we have the initial condition in M as 

d a M(x,0 + ) 


dx* 


- -f(x) 


( 21 ) 


For further simplicity, we also assume that the beam is simply supported at 
both ends and that the lateral load vanishes at both ends. Due to the 
symmetric nature of the problem as previously mentioned, the boundary 
conditions follow as 

£M e»M 

dx = ST** = 0 at *=0 (22a) 


d*M 

^7 = M=0 at x=L (22b) 


i 


It is important to point out that the neutral axis does not coincide with 
the centroidal axis in this beam problem since the viscosity is inhomogeneous 
due to its dependence on a non-uniform temperature distribution [36]. Since 
the axial force is zero in this problem, the distance to the neutral axis e Q 
may be determined in the first stage of damage from 



(23) 


It is our task now to extend the above mathematical formulation, which is 
valid only for the first stage of creep damage, into the second stage of creep 
damage. The shear stresses in technical beam theory are usually negligibly 
small when compared with the axial stress. It is thus reasonable to utilize 
the uniaxial strain -controlled damage law. The creep damage then follows from 
eqn. (12), which with the use of eqn. (17) yields 


t 

D(x.z.t) = c 0 jj [— ] n U-e 0 )dfj (24) 
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Experimental evidence [37] has shown that there is virtually no creep damage 
in a crystalline material under compression. Therefore, the above equation is 
valid only in the region e Q <z<h Q (see Fig. 1c), while the creep damage is 
assumed to be identically zero in the remainder of the region. The initial 
rupture time t^ may now be obtained from the implicit relation 


D 


cr 



(25) 


where the initial rupture clearly occurs at the midpoint of the bottom fiber 
(i.e. x=0 and z=h ), since it is there that the tensile strain is maximum in 
magnitude. 

Rupture thus starts at the point x,z=0,h o , and then develops into a 
moving front which in turn causes the beam to thin (see Fig 1 c ) . We shall 
call the region 0<x<6(t) the thinning zone, and the remaining interval 
6 ( t: ) <x<L the uniform zone since this interval is of uniform thickness. The 
quantities h and e, which designate the thickness and the distance to the 
neutral axis within the thinning zone of the beam, are clearly function of x 
and t. Furthermore, the flexural rigidity $ in the thinning zone is also a 
function of x and t since it involves h and e. Governing equation (18) with P 
given by eqn. (20) may now be restated in the thinning zone as 

^ - 0 0<x<Mt>. tj<t <26a) 

and in the uniform zone as 

+ k( ?L ) 1 = o 6(t)<x<L, tj<t (26b) 

dz 4 dt $ Q 
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where 


Mx.t) 1 

jd-JlKx.t) -b J [ Z *'dz'» 0£x<8(t). tjit (27) 

o 

It Is readily seen that governing equations (26) are subjected to a 
moving junction, which separates the thinning zone from the uniform zone. 
Note that the upper limit h(x,t) and the quantity e(x,t) In Integral (27) are 
changing and unknown functions, and thus we must obtain conditions which 
govern the variables h,4, and e. It may be shown that If the rupture front £ 
is prescribed as 

E : z = Mx.t) tjit 

eqn. (15) can be rewritten as 

dD dh 3D 

H * « H ■ 0 tin (28) 

Substitution of the expression for damage [eqn. (24)] into the above equation 
yields after some manipulation 

t 

^ = -(J) n (b-e) / | J (|) n d t' | . 0<x<6(t). tj<t (29) 

o 

The creep damage at the junction point Q in Fig. lc with coordinates 
x=6(t) and z=h Q should be equal to the critical value, i.e. 

D(x.z.t) * ® cr * x=6(t), z=h Q , tj<t 
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Following the same procedure employed for eqn. (29), the total time derivative 
of this equation now gives 


dMt) 

dt 


-M n / 



x*6(t) , 


tji.t 


(30) 


Note that the quantity «^(x,t) does not appear in the above equation as It Is a 
constant at the junction point Q (see Fig. 1c). 

Finally, eqn. (23) for the distance to the neutral axis of the beam now 
becomes for the second stage of creep damage 

h(x.t) 1 

f r z*-e(x, t) -i n 

J \- ~Tu 7 ) J dz ' " 0 ' 0<x<5(t). t x <t 


Differentiating the above equation with respect to time we obtain the equation 

I . 
n 

(z*-c) 

i 

o n 

[A( z ' ) 1 

where A(h) is the reciprocal viscosity function A(z) evaluated at z=h. 

We have thus obtained governing equations (26) subjected to Interface 
equations (29), (30) and (31), and we must solve these equations for the 

unknowns M,h,4, and e with boundary conditions (22). Although boundary 

conditions (22) are not applied at a moving boundary, we do have the interface 
equations which are applied at the moving junction. We shall thus use 

familiar terminology and call this problem a Stefan-like problem. It is 

readily seen that the present nonlinear problem is very complicated and 
numerically challenging; we shall present its numerical solution in the next 
section. 


dz'j. Oi: 


x<6(t) 


ti<t 


(31) 
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4. 


NUMERICAL TECHNIQUE ANO RESULTS 


4.1 Solution Technique 


For convenience we introduce the following nondimensional variables: 


' ■ T- • 


e = 


5 “ E 


x = £ . (Oixil) 


D = - — . (0<B<1) 

D cr 


z = ~ , (0<I<1) 

h o 


X = 


A( z) 
A( z=0) 


(32) 


E - ^ • <E 0 =1> 


•9 


_ 


M 


P o L * 


In the above, represents the flexural rigidity of a beam at a uniform 
perature T , where T y is 
present non-uniform beam. Thus 


temperature T , where T y is the temperature at the upper surface of the 


n +2 


/- 


nbh 


2(2n+l) [2A( z=0) ] 


(33a) 


and 




(l+2n)2 


n +1 h 


n z 9 -t -in _ _ 

— — ] z'dz' 

« 5(2') 


(33b) 


We now follow the practice that unless otherwise noted all variables without 
bars appearing in this section from this point on will be dimensionless 
variables. In accordance with the above definitions, governing equation (18) 
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may be reformulated in terms of dimensionless variables as 


a»M 

dx«at 


+ B(*f) n = o 


(34) 


In which we have Introduced the dimensionless quantity 


B 


Kt I P|J - 1 L 2n+2 


(35) 


The variables appearing on the R.H.S. of eqn. (35) are all In dimensional 
_ * 

form, and^ o Is evaluated by setting h=h Q In eqn. (33a). Note that the 

quantity B will be the key parameter In the present nondimenslonal study. 

Employing the same techniques presented In [35], we may eliminate the 
spatial partial derivative appearing in eqn. (34) for the first stage. We 
thus obtain the Integro-differentlal equations 



(36a) 


where 


F(x.x') - -<x-x')* (36b) 

G(x.x') = 3x* - 3xV - x'* + 3x'* - 2 ( 36c ) 

The numerical solution to the above equation with geophysical data for a beam 
of constant thickness was presented In [35]. Here, we will employ a numerical 
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technique other than the one presented in [35]. A discretization scheme using 
the method of lines in space is obtained from the above eqn. (36a). It 
follows that 



(37a) 


(37b) 


where 

x £ = (i-l)Ax = (i-l) 

n 2 

Note that the first integral in eqn. (37a) vanishes at 1=1, eqn. (37b) 
corresponds to boundary condition (22b), and N ? designates the number of 
spatial increments. Evaluating the Integrals by the Newton-Cotes formulas, we 
thus obtain a system of OOE's which may be solved by Gear's stiff ODE 
algorithm [38]. The result obtained above furnishes the solution in the first 
stage of damage, and provides the Initial data for the second stage of damage. 

Returning now to eqns. (26), we may again Integrate out the spatial 

derivatives to obtain for the thinning zone in the second stage 

x 8<t) 

£-?{/ F<x.x-)(j£>»dx' + J 0(,.x’)(|)"dx- 
o o 

1 M 1 

+ J G(x.x') (— ) n dx' | , 0<x<6( t) , lit (38a) 

lTa ' 
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and for the uniform zone 



where $ was defined in (33b), and F(x,x'), G(x,x') were defined in eqns. (36 
b,c). Note that we have used continuity of M and its derivatives at the 
junction point, e.g. M(x=4(t) + ,t) = M(x=6(t) ,t). In other words, the bending 
moment, shear force, deflection, and slope of the beam are all continuous at 
the junction point. 

In order to mathematically fix the moving junction and the limits of 
integration appearing in eqns. (38), we employ the concept of Landau's 
transformation [39] and introduce the variable changes 


X 

S = . ■■ , for thinning zone 0<x<6(t) 

6( t ) 


(39a) 


x-6(t) 

n = 7ZTT7\ * for mi f orm 6(t)<.xil (39b) 


Under such a transformation the partial time derivative a( )/at is replaced 
by the substantial time derivative D( )/Dt in accordance with 

8( ) D( ) dx 9( ) 
dt Dt dt dx 
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With the use of the chain rule, the transformed governing equations are 
obtained for the thinning zone as 


DM r _ d8(t) 3M g f 

Dt = 6(t) dt 71 Z l 


- + ? {set) J F( ? . ? '>(J> n j GU.C')<J> n dC' 


1 M 1 

+[ 1—8 ( t ) 3 J G( ? .n')(jj-) n dq' | . 0<C<1. ut 


(40a) 


and for the uniform zone as 

DM l-n dMt) 3M fi / \ M „ l M 

* ■ — a + I l 6< *> J J F(,.nO( r )"d V 


1 1 
+6(t) J G(Ti, c ')(J) n dC'+[l-6(t)] | G(n.V)(^-) n dT,'| . 


0<T1<1. 


(40b) 


Similarly, an application of Landau's transformation to interface 
equations (29), (31) yields 


Dh _ z d6(t) 3h 
Dt = 8(t) dt 


(J) n (h-e) / | J[«] n dt‘| . 0 <?< 1 . l<t 


t; ^8(t) de rDh ^ d6(t) dh-irh— c *jn 

&(t) dt 3 c +n *-Dt 8(t) ~dt il-ll-A(li) J / 

1-1 

r r (z ' _e) i 

| J F dz J ' 01 * llt; 

o n 

[A(z') J 
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The transformed interface equation (30) has the slightly different form 


d6(t) 

dt 


= -6<t)M n / jn J M n_1 ™ dt') • C =1 • U* 


u b | 

ac J 


(43) 


since 4(t) involves only the single variable t. By virtue of variable 

changes (39), the substantial (material) time derivative D( ) /Dt appearing in 

the above transformed equations possesses a numerical value identical with 

[3( )/at]„ or [a( )/3tJ . It should be noted that fixing the moving junction 
C n 

unfortunately leads to governing equations of even more complicated form, and 
the above set of equations surely will not have a closed form solution. A 
numerical scheme will now be presented. 

First, the method of lines in space [40] is utilized to eliminate the 
spatial derivatives from the above integro-dif ferential equations in 
accordance with the discretization scheme shown in Fig. 2. Accordingly, we 
employ the central finite difference approximations for the interior points 
and one-sided three point formulas for the end points A,B and the junction 
point Q (see Fig. 2). Accordingly, eqns. (40-43) yield the followina: 


dS( t) 
dt 


= -2AC 6(t) M? / jn J Mf' 1 tM._ 2 -4M._ 1 +M i ]dt' J . 


i=N : +l (44) 


DM, 

Dt 


1 1 

i = | |S( t) J G(q. £') (J) n dc'+[l-6(t)l j G( q.n') (^-) n dTi* j, i=l (45a) 
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DM, 


’i d 6 ( t) 


Dt 2 A 5 6 ( t ) dt 


[M i . 2 - 4 M i _ 1 + 3 M i ]+ | (s<t) J F( Ci . ? ')(|) n d C ' 


1 1 Ml 

+ 6 (t) J GU 1# ;’)(j*) n di;'+[l- 6 <t)] J GtCj.n') (^-) n dtj'J . i-Nj +1 


5T - g ^mi TT 1 ‘WW- ? {“*> J *<m- ?'> 


1 

+ [ 1 - 8 ( t ) ] J F(t li .ti') (J-) n dTi'+ 6 (t) J G( ni ,C') (|) n dc* 

0 ° o 

1 

+ [ 1 - 8 ( t ) ] J G( ni .^') (j-) n di,' j. i=N 1 + 2 ,N 1 + 3 ,....N 1 +N 2 


DM, 


Dt 0 ' i-Nj+Nj+l 


Dh. M, t - m. 

— , J ; . 1 , I 

o 


M. 


Dt ^ ' 1 ' l J 


Dhj 

Dt~ 


C i d 6 (t) M i f \ M i 'i 

2 A? 6 ( t ) “dt - [h i+l" h i-l 1 -( ^ )n(ll r e i ) 7 { J ( t p7 )n<lt, J ’ 


i= 2 ,3 . . . . ,N, 


Dh j 

Dt 


= 0 , 


i=N 1+ l 


( 45 C) 


( 45 d) 


( 45 e) 


( 46 a) 


( 46 b) 


( 46 c) 
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t - 1 


Dei D hi hi-ti 1 > i (z'-ej) 

Dt n Dt L A (ji i )J ^ \J 2. 


dz 1 


i=l 


(47a) 


[A(z') ] 


Dei 


? i d6( t ) 


Dt 2AC 6(t) dt 


r V e i,n f f 1 1 

/{/ r** 1 )- “ 2 - 3 »i < 4,| » 


Dhi 

Ie i+l“ e i-l ] + 4^- 


^ d6(t) (h _ h vl 

2AC 6(t) dt (hi+1 J 


1-1 


o n 

[A(z' ) ] 


D#i 

Dt 


= 0 . 


l=N 1+ l 


(47c) 


In the above 


- (i-l)AC = — (i-1) , i=l,2 .... ,N, 


r\i = (i-Nj-DAr, = — ( i-Nj-1) , 
N 2 




i=Nj+l, N x +2 Ni+N 2 +1 


and N 1 ,N 2 designates, respectively, the number of spatial Increment in 
thinning zone and uniform zone, and Jf. is obtained by setting h=h^ is eqn. 
(33b). 

Note that eqn. (45e) corresponds to boundary condition (22b). Moreover, 
Dh/Ot and De/Dt as given by eqns. (46c) and (47c) vanish at the junction 
point, since at any instant we always have h=l and e=e Q at this point. The 
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integrals appearing in the above set of equations were evaluated by use of the 
Newton-Cotes formulas. We thereby obtained a large system of 3N 1 +N 2 +4 
ordinary differential equations. A computer program was developed first to 
solve eqn. (23) for e Q and then to solve eqns. (37) for the M^'s using the 

initial bending moment function (see Sec. 4.2) as the initial condition. We 
then used these results along with h..=l, e.=e Q , 4=0 as the initial conditions 
to solve the system of OOE's, eqns. (44-47). 

It is also useful to compute the deflection of the beam. This can be 

accomplished by substituting the above bending moments M. into the 
nondimensinal form of the equation of equilibrium 
2 

a M 

= -P + Kw 

7 ~ 

3 X 

However, this approach leads to considerable numerical error due to the 

presence of the derivative term in the above equation. An alternate and more 
accurate approach is to develop differential equations in the deflection. The 
same solution technique used for the bending moment is also applicable to 

these differential equations in the deflection. For the sake of brevity, the 
resulting equations will not be presented here (see [18]). An even larger 
system of 4N 1 +2N 2 +5 ODE's is obtained in this case. A high accuracy yet 
costly numerical algorithm, i.e. Gear's stiff 00E algorithm, was used to solve 
this system. 
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4.2 Solutions and Discussion 

Our attention Is first directed to a special case In which closed form 
solutions exist. Thus, let us delete the elastic Winkler foundation and also 
consider a beam with a uniform temperature distribution equal to T y (a 
dimensional quantity). Under such circumstances, the bending moment M remains 
constant in time, and the neutral axis coincides with the centroldal axis 
owing to the homogeneous nature of the material properties. Moreover, we have 
the following values for the dimensionless quantities [see eqns. (32) and (33)] 


A(z) =1 

« +2 
$ = h“ 

The dimensionless governing equations for h,e [see eqns. (29) and (31)] in the 
thinning zone (0 <x<6(t)) follow for this special case as 



at t 

2 J h“ 1-2n dt' 


o 


Ut 


de „ dh 

at " 2 at 


lit 


(48a) 


(48b) 


and that for 4(t) becomes 


d6(t) 

dt 


M 



x*6 ( t ) , lit 


(48C) 
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Note that these equations may be solved consecutively and that eqns. (48 
a,b) include neither the x variable nor the input load function P(x,t) 
explicitly. Physically, eqn. (48b) indicates that although 2e 0 =h 0 =1 
initially, both quantities will be equal to zero at the instant the beam 
collapses. After eliminating the integral via differentiation, eqn. (48a) may 
be rewritten as 


h $ • » 2 - 2 ] - » < 49 » 

This differential equation may be solved analytically with the initial 
conditions 


h = l at t = tj 


and 


ah l 

Hi = ~ 2t^ 8t * = *1 

where tj=tj(x) designates the time required for a material point with 
coordinates (x,l) to reach the critical state. The second initial condition 
in the above was obtained by setting t=tj and h=l in eqn. (48a). The solution 
to eqn. (49) is then obtained as 

4 2 2 . 1+2 a 

• # ^<«(‘> • lit (50) 

which is identical to the result derived in [8]. 
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The solution to eqn. (48b) for e then follows directly from the solution 
of eqn. (50) for h. Here tj(x) may be expressed in terms of the bending 
moment M(x). Since the point (0,1) reaches the critical state at time t=l 
while the point (x,l) ruptures at time t=tj(x), it follows from eqn. (25) that 


M 

t x <x) = < r >* 


(51) 


Here, M Q and M represent the bending moments at points with x-coordinates 
equal to 0 and x, respectively. Equation (50) thus becomes 


(Vt = — h 2 ”- 1 - — . 0<x<Mt) . l<t <52) 

M e 1-2 n l-2n 

Although differential equation (48a) in h does not explicitly involve the 
variable x and the bending moment M, its solution (52) is seen to be directly 
related to M. 

The function t j ( x) in eqn. (51) may be inverted in the simple case that 
the bending moment M(x) is monotonic in nature. In fact, for such a simple 
case the constraint on the moving juntion 

x = 8(t) 

is invertible and is physically equivalent after inverstion to 

t - tj(x) 

Consequently, witht he use of eqn. (51), eqn. (48c) attains the alternate form 


d8(t) 

dt 


M n+1 


_un 25 

ax 


x=6( t) . l<t 


(53) 
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Consider for the moment a very simple case in which a point-load is applied on 
the beam at x=0. The bending moment for this load is simply given as 

M = M o (i- x ) (54) 

and differential equation (53) becomes 

^ - i ll-Ktl]** 1 . lit 

With the initial condition 4 at t=l equal to 0, the above equation yields the 
solution 

_ l 

6(t) - 1 - t * . l<t (55) 

Note that solution (55) is also obtainable from the solution (52) by setting 
h=l and using expression (54). The solutions of this special closed-form case 
are now complete. 

We thus turn our attention to the original problem, containing in general 
both the Winkler term and a temperature gradient. A singularity appears in 
eqn. (39a) at time t=l when 4=0, and this leads to numerical difficulties. 
This obstacle may be circumvented by Introducing an "imperfection" [33] in 
eqn. (39a). Here, we follow the latter approach and Introduce an imperfection 
in 4 as 

6(t) = 1.0x10'* at t=l 

Moreover, the temperature in the beam is assumed to be linearly distributed In 
the z-direction in accordance with (see Fig. lc) 

T " T u + (T b " V z * °^ 1 
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where the dimensional surface and bottom temperatures are chosen respectively 
as 

T n = 300°k . T b = 360°k 

Also, we use for the creep activation energy the dimensional value 
AH = 0 .112x10* J-mole -1 

Note that because of the nondimensional form of our governing equations it was 
not necessary to stipulate a specific material. Finally the number of 
increments chosen in the present beam problem were 

* 5 for the thinning zone 
N 2 = 10 for the uniform zone 

which yield a total of 29 ODE's, or 45 if the equations for deflection are 
also included. In order to limit the complexity of this nonlinear problem, 
only the uniformly applied load is considered here. In this case, the shape 
function [see eqn. (20)] of the applied load is simply 

f(x) • l , 0<x<l 

The initial bending moment function M(x,0 + ) is obtainable from eqns. (21) and 
(22) , and is given as 

M(x,0 + ) = (x*-l) / 2 (56) 

which will be used as the initial condition for eqns. (37). 

The results we shall present may be separated into two groups, i.e. those 
with the foundation present and those with the foundation absent. In the 
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latter case of the foundation absent, we have from eqn. (35) B=0 since in this 
case the spring constant for the foundation is identically zero. We consider 
the B=0 case first, and note that here the bending moment is independent of 
time, and is thus given simply by eqn. (56) after the lateral load is 
applied. (For the B=0 case, we did not calculate the deflection of the 
beam.) Figures 3 and 4 display the propagating rupture front for 8=0 in the 
second stage of damage for n=3 and 5, respectively. In these figures the 
depth and axial coordinates z,x of the beam are given in nondimensional form. 
The sequence of curves inside the beam trace the propagation of the rupture 
front as time ellapses. The S(t) function at time t is given by the distance 
along the bottom surface (z=l) from the point x=0 to the intersection of the 
curve for time t with the bottom surface. Note that the beam of n=3 material 
exhibits a wider thinning zone than does the beam of n=5 material. It would 
appear that the rupture front of the n=5 beam propagates faster than that of 
the n=3 beam. But you are reminded that this observation is made for a 
nondimensional time scale and will not necessarily follow for dimensional 
time. Each set of these curves for a parameter n required about 1 minute of 
computer time (CPU) . 

We now turn to the general case with the Winkler foundation present. The 
dimensionless parameter B [see eqn. (35)] contains a group of constants 
including the spring constant of the foundation, applied load, geometry and 
material properties of the beam, and is considered arbitrary in the present 
nondimensional study. Here we chose for illustration the value B=1 which 
requires that the numerator and denominator in eqn. (35) be of the same 
order. In addition to the bending moment M, we also computer the deflection w 
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for this case. Due to the complicated nature of the governing large system of 
ODE's, a considerable amount of computer time was required to complete one run 
for a specific value of n. Thus we limited the computer time to 1000 seconds 
(about 16.7 minutes) per run, and accordingly obtained a reasonable number of 
solution curves for time steps in the early part of the second stage of 
damage. The computations could have easily been extended up to the point of 
final collapse of the beam, but for reasons of economy this was not done. 

Figures 5a, b give respectively the nondimens ional deflection of the beam 
for values 3 and 5 for the stress power n. Since the chosen load is uniformly 
distributed, these curves do not exhibit the characteristic "uplift" which 
often occurs under centrally concentrated loads or point loads on a beam with 
a Winkler foundation [34, 35]. According to the flexural model presented in 
[35], the deflections of a beam which experiences no damage approaches an 
asymptotic limit as the time tends to infinity. Howler, no asymptotic 
deflection solution exists in the present problem, since damage causes the 
beam to thin and accordingly the deflection is unbounded. This may readily be 
seen in Fig. 5b, in which the increment of beam deflection is clearly 
increasing in the final few time steps shown. Although we have used Norton's 
steady creep law, eqn. (16), to describe the creep behavior of the material, 
the nature of the deflection shown in Fig. 5b is similar in form to the 
typical creep curve with its three stages of creep. Such behavior coincides 
with the recent experimental investigations [41, 42] in which a beam with a 
deep notch was subjected to a uniform temperature and point load. This can be 
explained by the fact that as the beam starts thinning the remaining material 
carries the same load but with greater stress. 
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In [35] where no damage was included, McMullen et al noted that the 
bending moment relaxes after the load is applied and approaches zero as time 
tends to infinity. Figures 6a, b exhibit this same relaxation of the bending 
moment in the more general case where damage causes the beam to thin. 
Furthermore, Fig. 6b shows that the relaxation of the bending moment 
accelerates in the final few time steps shown; it is believed that Fig. 6a 
would also do the same if the computer time had been extended. Since the 
lifetime of the beam is finite, the beam should collapse before the bending 
moment vanishes. Figures 7a, b display the propagating rupture front for 8=1 
with n=3 and 5, respectively. As in the B=0 case, the rupture fronts in the 
present case have sharper profiles in the n=5 beam than in the n=3 beam. And 
the rupture front for the n=5 beam propagates faster than that for the n=3 
beam relative to the nondimensional time scale. We also point out that the 
numerical scheme for the system of ODE's presented in the previous section is 
stiffer for n=5 than that for n=3, since within the chosen limit of compute 
time (100 seconds) the final time step reached was t=1.90 for the case of n=3 
and only t=1.60 for the case of n=5. 

It should be noted that the results displayed may contain some numerical 
error in the later time intervals, since we are restricted by the limitations 
of infinitesimal strain and small rotations. Although we have formulated the 
problem in an idealized manner, a significant amount of mathematical 
difficulty was still encountered. If one attemps to relax some of the 
assumptions imposed, the complexity of the problem could increase greatly and 
possibly preclude a successful numerical analysis. 
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P(x,t) 




^ High temperature 
^ foundation 


Fig. la Beam on high temperature foundation, subjected to lateral 
load P(x,t). 



Fig. 1b Beam with simple end supports, elastic Winkler foundation, 
and symmetric load. 


177 



178 


Fig. 1c Propagation of the rupture front and the linear temperature 
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Fig. 4 Nondimensional deflection of a beam resting on Winkler 

foundation — (a) B=1, n=3, T b =360°k and T u =300°k, (b) B=1, 
n=5, T b =360°k and T u =300°k. 
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Fig 4. (continued) 
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(a) 

Fig. 5 Nondimensional bending moment along a beam on Winkler 
foundation — (a) B=1.0, n=3, T b =360°k and T u =300°k, (b) 
B=1.0, n=5, T b =360°k and T u =300°k. 
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Fig. 5 (continued) 



(b) 
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Fig. 6 Propagation of rupture front in a beam on Winkler foundation — 
(a) B=1, n=3, T b =360°k and T u =300°k, (b) B=1, n=5, T b =360°k 
and T u =300°k. 
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